Decomposing time series with complex seasonality

Rob J Hyndman

COMPSTAT: 23 August 2022

Complex seasonality

Complex seasonality

Complex seasonality

Complex seasonal topology

Example: hourly data

Complex seasonality

  • Multiple seasonal periods, not necessarily nested
  • Non-integer seasonality
  • Irregular seasonal topography
  • Seasonality that depends on covariates
  • Complex seasonal topology

No existing decomposition method handles all of these.

Two solutions

  1. MSTL: For multiple integer seasonal periods.
  2. STR: For all types of complex seasonality.

MSTL

 

  • Kasun Bandara, Rob J Hyndman, Christoph Bergmeir (2022) MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns. International J Operational Research, to appear. robjhyndman.com/publications/mstl/

  • Implemented in R packages forecast and fable.

MSTL

velec |>
  model(STL(Demand)) |>
  components() |>
  autoplot()

MSTL

MSTL

y_t = T_t + \sum_{i=1}^I S_t^{(i)} + R_t

y_t= observation at time t
T_t= smooth trend component
S_t^{(i)}= seasonal component i
i = 1,\dots,I
R_t= remainder component

Estimation

Components updated iteratively.

MSTL

# X: time series object
# periods: vector of seasonal periods in increasing order
# s.window: seasonal window values
# iterate: number of  STL iterations

seasonality <- matrix(0, nrow = nrow(X), ncol = length(periods))
deseas <- X
for (j in 1:iterate) {
  for (i in 1:length(periods)) {
    deseas <- deseas + seasonality[, i]
    fit <- model(
        STL(deseas ~ season(period = periods[i], window = s.window[i]))
      ) |>
      components()
    seasonality[, i] <- fit$season
    deseas <- deseas - seasonality[, i]
  }
}
trend <- fit$trend
remainder <- deseas - trend
return(trend, remainder, seasonality)

MSTL

fable syntax

tsibble |>
  model(STL(variable) ~ season(period = a, window = b) +
                        season(period = c, window = d))



forecast syntax

vector |>
  msts(seasonal.periods = c(a, c)) |>
  mstl(s.window = c(b, d))

STR

 

  • Alex Dokumentov and Rob J Hyndman (2022) STR: Seasonal-Trend decomposition using Regression. INFORMS Journal on Data Science, to appear. robjhyndman.com/publications/str/

  • Implemented in R package stR.

STR

y_{t} = T_{t} + \sum_{i=1}^{I} S^{(i)}_{t} + \sum_{p=1}^P \phi_{p,t} z_{t,p} + R_{t}

T_t= smooth trend component
S_t^{(i)}= seasonal component i (possibly complex topology)
z_{p,t}= covariate with coefficient \phi_{p,t} (possibly time-varying)
R_t= remainder component

Estimation

Components estimated using penalized MLE

Smoothness via difference operators

Smooth trend obtained by requiring \Delta_2 T_t \sim \text{NID}(0,\sigma_L^2)

  • \Delta_2 = (1-B)^2 where B= backshift operator
  • \sigma_L controls smoothness

f(\bm{D}_\ell \bm{\ell}) \propto \exp\left\{-\frac{1}{2}\big\|\bm{D}_\ell \bm{\ell} / \sigma_L\big\|_{L_2}^2\right\}

  • \bm{\ell} = \langle T_{t} \rangle_{t=1}^{n}
  • \bm{D}_\ell= 2nd difference operator matrix: \bm{D}_\ell\bm{\ell} = \langle\Delta^2 T_{t}\rangle_{t=3}^n

Smooth 2D seasonal surfaces

  • m_i= number of “seasons” in S^{(i)}_{t}.
  • S^{(i)}_{k,t}= 2d season (k=1,\dots,m_i;t=1,\dots,n)
  • \sum\limits_k S^{(i)}_{k,t} = 0 for each t.

Smooth 2D seasonal surfaces

  • \bm{S}^{(i)} = [S_{k,t}^{(i)}] the ith seasonal surface matrix
  • \bm{s}_i = \text{vec}(\bm{S}_i)= the ith seasonal surface in vector form

Smoothness in time t direction:

\begin{align*} \bm{D}_{tt,i} \bm{s}_i &= \langle \Delta^2_{t} \bm{S}^{(i)}_{k,t} \rangle \sim \text{NID}(\bm{0},\sigma_{i}^2 \bm{\Sigma}_{i})\\ f(\bm{s}_i) &\propto \exp\Big\{-\frac{1}{2}\big\|\ \bm{D}_{tt,i}\bm{s}_i / \sigma_i\big\|_{L_2}^2\Big\} \end{align*}

Analogous difference matrices \bm{D}_{kk,i} and \bm{D}_{kt,i} ensure smoothness in season and time-season directions.

Gaussian remainders

  • R_{t} \sim \text{NID}(0,\sigma_R^2).
  • \bm{y} = [y_1,\dots,y_n]'= vector of observations
  • \bm{Z}=[z_{t,p}]= covariate matrix with coefficient \bm{\Phi} = [\phi_{p,t}]
  • \bm{Q}_i= matrix that extracts \langle S^{(i)}_{\kappa(t),t} \rangle_{t=1}^{n} from \bm{s}_i.
  • Residuals: \bm{r} = \bm{y} - \sum_i\bm{Q}_i\bm{s}_i -\bm{\ell} - \bm{Z}\bm{\Phi} have density

f(\bm{r}) \propto \exp\Big\{-\frac{1}{2}\big\|\bm{r}/\sigma_R\big\|_{L_2}^2\Big\},

MLE for STR

Minimize wrt \bm{\Phi}, \bm{\ell} and \bm{s}_i:

\begin{align*} -\log \mathcal{L} &= \frac{1}{2\sigma_R} \Bigg\{ \Big\| \bm{y}- \sum_{i=1}^I \bm{Q}_i\bm{s}_i - \bm{\ell} - \bm{Z}\bm{\Phi} \Big\|_{L_2}^2 + \lambda_\ell\Big\|\bm{D}_\ell \bm{\ell}\Big\|_{L_2}^2 \\ & \hspace*{1cm} + \sum_{i=1}^{I}\left( \left\|\lambda_{tt,i} \bm{D}_{tt,i} \bm{s}_i \right\|_{L_2}^2 + \left\|\lambda_{st,i} \bm{D}_{st,i} \bm{s}_i \right\|_{L_2}^2 + \left\|\lambda_{ss,i} \bm{D}_{ss,i} \bm{s}_i \right\|_{L_2}^2 \right) \Bigg\} \end{align*}

Equivalent to linear model

\bm{y}_{+} = \bm{X}\bm{\beta} + \bm{\varepsilon}

  • \bm{y}_{+} = [\bm{y}',~ \bm{0}']'
  • \bm{\varepsilon} \sim N(\bm{0},\sigma_R^2\bm{I})

\bm{X} = \begin{bmatrix} \bm{Q}_1 & \dots & \bm{Q}_I & \bm{I}_n & \bm{Z} \\ \lambda_{tt,1} \bm{D}_{tt,1} & \dots & 0 & 0 & 0 \\ \lambda_{st,1} \bm{D}_{st,1} & \dots & 0 & 0 & 0 \\ \lambda_{ss,1} \bm{D}_{ss,1} & \dots & 0 & 0 & 0 \\ 0 & \ddots & 0 & 0 & 0 \\ 0 & \dots & \lambda_{tt,I} \bm{D}_{tt,I} & 0 & 0 \\ 0 & \dots & \lambda_{st,I} \bm{D}_{st,I} & 0 & 0 \\ 0 & \dots & \lambda_{ss,I} \bm{D}_{ss,I} & 0 & 0 \\ 0 & \dots & 0 & \lambda_\ell \bm{D}_{tt} & 0 \end{bmatrix}

STR

Three seasonal components, quadratic temperature regressors

STR outliers


For more information

Slides: robjhyndman.com/seminars/compstat2022

Find me at:

robjhyndman.com
@robjhyndman
@robjhyndman
rob.hyndman@monash.edu